% Calculate Fermi level of intrinsic InSb for parabolic bands.
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM.
% Based on code written by Carlos A. Armenta and Sonam Yadav.
% Date: 18 October 2023
% Required package: Symbolic Math Toolbox Add-On. 
%
% References:
% Melissa Rivero Arias, Carlos A. Armenta, Carola Emminger, Cesy M. Zamarripa, 
% Nuwanjula S. Samarasingha, Jaden R. Love, Sonam Yadav, and Stefan Zollner
% J. Vac. Sci. Technol. B 41, 022203 (2023).
%
% Physical constants (energy units in eV unless noted) from CODATA at NIST
hbar=6.582119569E-16; % eV*s
e=1.602176634E-19; % A*s
m0=9.1093837015E-31; % kg
eps0=8.8541878128E-12; % F/m
kB=0.08617333262E-3; % eV/K Boltzmann constant
kBJoule=kB*e;
c=299792458; % m/s speed of light in vacuum
%
h=hbar*2*pi; % eV*s
hJoule=h*e; % J*s
hbarJoule=hbar*e;
hbarc=hbar*c*100; % eV*cm
hbar2m0=hbar*hbarJoule/m0*1E20; % hbar^2/m0 in Angstrom^2
Rydberg=e^4*m0/(8*eps0^2*hJoule^2)/e; % eV, hydrogen ground state energy
Hartree=2*Rydberg; % eV, atomic unit for energy 
aBohr=(hJoule^2*eps0)/(pi*e^2*m0)*1E10; % Angstrom
%Prefactor e^2/(eps0*m0) for plasma frequency
PlasmaFactor=e^2/(eps0*m0)*1E6; % cm-3/s^2
DebyeFactor=100*e^2/(eps0*kBJoule); % Kcm
% amplitude prefactor for optical interband transitions
Afactor=4*sqrt(Rydberg)/3;
% prefactor for the electron and hole concentrations in cm^{-3}
Nfactor=((2*m0*kB)/(pi*hbar*hbarJoule))^1.5*1E-6/4;
%
% Material parameters for InSb (compare Ioffe Institute web site)
a=6.47931; % lattice constant of InSb in Angstrom, from Adachi (Springer)
% E0=0.237; % eV, from Lawaetz 1971
E0=0.2; % eV from supplemental materials in Melissa's paper
Delta0=0.81; % eV, from Lawaetz 1971
EP=23.1; % eV, momentum matrix energy, from Lawaetz 1971
% EP=0.44*Hartree*2; % convert 2P^2 from atomic units to eV, from Kane 1957
% me=1/(1+EP/3*(2/E0+1/(E0+Delta0))) % should be about 0.014 (Ioffe)
me=0.014; % froim Ioffe Institute web site
% mlh=3*E0/(2*EP-3*E0) % should be about 0.015 (Ioffe)
mlh=0.015; % from Lawaetz 1971
mso=1/(EP/3/(E0+Delta0)-1); % should be about 0.19 (Ioffe) or 0.107 (Lawaetz 1971)
% mhh=0.47; % from Lawaetz 1971
mhh=0.43; % from Ioffe Institute web site 
% mso=0.107; % from Lawaetz 1971
eps_infty=15.7; % high-frequency dielectric constant (Ioffe)
%Calculated parameters
muhh=(me*mhh)/(me+mhh);
mulh=(me*mlh)/(me+mlh);
Rhh=muhh/eps_infty^2*Rydberg; % eV
Rlh=mulh/eps_infty^2*Rydberg; % eV
ahh=eps_infty/muhh*aBohr; % Bohr radius in Angstrom
alh=eps_infty/mulh*aBohr; % Bohr radius in Angstrom
%
% Calculate Fermi energy for parabolic bands. This is an approximation. 
% The conduction band and light hole band of InSb are not parabolic. 
%
T=300 % K (temperature)
Fermi=@(m,x)-polylog(m+1,-exp(x)); % Fermi-Dirac integrals expressed as polylogs.
phh=@(x)Nfactor*(mhh*T)^1.5*Fermi(0.5,-x/(kB*T));
plh=@(x)Nfactor*(mlh*T)^1.5*Fermi(0.5,-x/(kB*T));
n=@(x)Nfactor*(me*T)^1.5*Fermi(0.5,(x-E0)/(kB*T));
p=@(x)phh(x)+plh(x);
%
muvar=0.1:0.001:0.3;
% This plot shows the electron and hole densities versus chemical potential
% for a given temperature defined above. Compare Fig. S5 in Rivero Arias (2023). 
f=figure('Name','Electron and hole densities versus chemical potential', ...
    'Color','white');
f.Position(3:4)=[700 420];
pl=plot(muvar,n(muvar)/1E15,'r-',muvar,p(muvar)/1E15,'b-');
xlim([0.14 0.18]);
ylim([4 32]);
ax=gca; ax.FontSize=12; ax.YTick=4:2:32; ax.XTick=0.14:0.01:0.18; 
xlabel('chemical potential \mu (eV)'); ylabel('carrier density (10^{15} cm^{-3})'); 
legend('electron density','hole density','Location','northeast'); legend('boxoff');
pl(1).LineWidth=2; pl(2).LineWidth=2;
text(0.05,0.6,'T='+string(T)+' K','FontSize',14,'Units','normalized');
%
diff=@(x)n(x)-p(x);
muroot=fzero(diff,0);
murootmeV=muroot*1000
proot=p(muroot)
nroot=n(muroot)
text(0.05,0.5,'E_F='+string(round(murootmeV,3,"significant"))+' meV','FontSize',14,'Units','normalized');
%
% Plot the same values with a logarithmic vertical axis, similar to Fig. S6 in Rivero Arias (2023). 
f=figure('Name','Electron and hole densities versus chemical potential', ...
    'Color','white');
f.Position(3:4)=[700 420];
pl=plot(muvar,n(muvar)/1E15,'r-',muvar,p(muvar)/1E15,'b-');
xlim([0.14 0.18]);
ylim([4 34]);
ax=gca; ax.FontSize=12;
ax.YScale="log"; ax.YTick=5:5:32; ax.XTick=0.14:0.01:0.18; 
xlabel('chemical potential \mu (eV)'); ylabel('carrier density (10^{15} cm^{-3})'); 
legend('electron density','hole density','Location','northeast'); legend('boxoff');
pl(1).LineWidth=2; pl(2).LineWidth=2;
text(0.05,0.7,'T='+string(T)+' K','FontSize',14,'Units','normalized');
text(0.05,0.6,'E_F='+string(round(murootmeV,3,"significant"))+' meV', ...
    'FontSize',14,'Units','normalized');
%
% return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now calculate the Fermi energy as a function of temperature from 30 to 800 K.
% Compare Fig. S7 in Rivero Arias (2023). This is for parabolic bands.
Tstart=30;
Tstep=10;
for i=1:1:78
    Ta(i)=Tstart+(i-1)*Tstep;
    T=Ta(i);
    phh=@(x)Nfactor*(mhh*T)^1.5*Fermi(0.5,-x/(kB*T));
    plh=@(x)Nfactor*(mlh*T)^1.5*Fermi(0.5,-x/(kB*T));
    n=@(x)Nfactor*(me*T)^1.5*Fermi(0.5,(x-E0)/(kB*T));
    p=@(x)phh(x)+plh(x);
    diff=@(x)n(x)-p(x);
    muroot=fzero(diff,0);
    mupara(i)=muroot; % chemical potential (degenerate) with parabolic bands
    munondeg(i)=E0/2+3*kB*T/4*log(mhh/me); % chemical potential (non-degenerate)
    npara(i)=n(muroot); % carrier concentration (degenerate) with parabolic bands
    nnondeg(i)=Nfactor*(me*mhh)^0.75*T^1.5*exp(-E0/(2*kB*T)); % nondeg carrier concentration
    fprintf('T=%5.0f, E_F=%5.3f eV, n=%8.1e,  n=%8.1e\n',Ta(i),mupara(i),npara(i),nnondeg(i));
end
%
f=figure('Name','Chemical potential versus temperature (parabolic)', ...
    'Color','white');
f.Position(3:4)=[700 420];
pl=plot(Ta,mupara,'k-',Ta,munondeg,'k--');
xlim([0 800]); ylim([0.10 0.30]);
legend('degenerate parabolic','non-degenerate parabolic','Location','northwest'); 
legend('boxoff');
ax=gca; ax.FontSize=12; ax.YTick=0.10:0.02:0.30;
ylabel('chemical potential \mu (eV)'); xlabel('Temperature T (K)'); 
%
% Draw carrier concentration versus temperature with logarithmic vertical axis.
% Compare Fig. S8 in Rivero Arias (2023). This is for parabolic bands.
f=figure('Name','Carrier concentration versus temperature (parabolic)', ...
    'Color','white');
f.Position(3:4)=[700 420];
pl=plot(Ta,npara,'r-',Ta,nnondeg,'k--');
xlim([0 800]);
ylim([1E13 1E18]);
ax=gca; ax.FontSize=12; ax.YScale="log"; 
ylabel('carrier concentration n (cm^{-3})'); xlabel('Temperature T (K)');
legend('degenerate parabolic','non-degenerate parabolic','Location','northwest'); 
legend('boxoff');
pl(1).LineWidth=2; pl(2).LineWidth=2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
